Go back to the Contents page.
Press Show to reveal the code chunks.
# Set seed for reproducibility
set.seed(1982)
# Set global options for all code chunks
knitr::opts_chunk$set(
# Disable messages printed by R code chunks
message = FALSE,
# Disable warnings printed by R code chunks
warning = FALSE,
# Show R code within code chunks in output
echo = TRUE,
# Include both R code and its results in output
include = TRUE,
# Evaluate R code chunks
eval = TRUE,
# Enable caching of R code chunks for faster rendering
cache = FALSE,
# Align figures in the center of the output
fig.align = "center",
# Enable retina display for high-resolution figures
retina = 2,
# Show errors in the output instead of stopping rendering
error = TRUE,
# Do not collapse code and output into a single block
collapse = FALSE
)
# Start the figure counter
fig_count <- 0
# Define the captioner function
captioner <- function(caption) {
fig_count <<- fig_count + 1
paste0("Figure ", fig_count, ": ", caption)
}capture.output(
knitr::purl(here::here("functionality.Rmd"), output = here::here("functionality.R")),
file = here::here("purl_log.txt")
)
source(here::here("functionality.R"))Let \(\Gamma = (\mathcal{V},\mathcal{E})\) be a metric graph. Let \(u_d: \Gamma \times(0, T) \rightarrow \mathbb{R}\) be the desired state and \(\mu>0\) a regularization parameter. We define the cost functional
\[\begin{equation} \label{eq:costfun} \tag{1} J(u, z)=\frac{1}{2} \int_0^T\left(\left\|u-u_d\right\|_{L_2(\Gamma)}^2+\mu\|z\|_{L_2(\Gamma)}^2\right) dt \end{equation}\]
Let \(f:\Gamma\times (0,T)\rightarrow\mathbb{R}\) and \(u_0: \Gamma \rightarrow \mathbb{R}\) be fixed functions. We will call them right-hand side and initial datum, respectively. Let \(\alpha\in(0,2]\) and \(z: \Gamma \times(0, T) \rightarrow \mathbb{R}\) denote the control variable. We shall be concerned with the following PDE-constrained optimization problem: Find
\[\begin{equation} \label{eq:min_pro} \tag{2} \min\; J(u, z) \end{equation}\] subject to the fractional diffusion equation \[\begin{equation} \label{eq:maineq} \tag{3} \left\{ \begin{aligned} \partial_t u(s,t) + (\kappa^2 - \Delta_\Gamma)^{\alpha/2} u(s,t) &= f(s,t)+z(s,t), && \quad (s,t) \in \Gamma \times (0, T), \\ u(s,0) &= u_0(s), && \quad s \in \Gamma, \end{aligned} \right. \end{equation}\] with \(u(\cdot,t)\) satisfying the Kirchhoff vertex conditions \[\begin{equation} \label{eq:Kcond} \tag{4} \mathcal{K} = \left\{\phi\in C(\Gamma)\;\middle|\; \forall v\in \mathcal{V}:\; \sum_{e\in\mathcal{E}_v}\partial_e \phi(v)=0 \right\} \end{equation}\] and the control constraints \[\begin{align} \label{control_constraints} \tag{5} a(s,t)\leq z(s,t)\leq b(s,t)\;\mathrm{a.e.}\;(s,t)\in\Gamma \times(0, T). \end{align}\]
The optimal variables \((\bar{u}, \bar{p}, \bar{z})\) satisfy
\[\begin{equation} \label{eq:maineqoptimal} \tag{6} \left\{ \begin{aligned} \partial_t \bar{u}(s,t) + (\kappa^2 - \Delta_\Gamma)^{\alpha/2} \bar{u}(s,t) &= f(s,t)+\bar{z}(s,t), && \quad (s,t) \in \Gamma \times (0, T), \\ \bar{u}(s,0) &= u_0(s), && \quad s \in \Gamma, \end{aligned} \right. \end{equation}\] and \[\begin{equation} \label{eq:adjointeq} \tag{7} \left\{ \begin{aligned} -\partial_t \bar{p}(s,t) + (\kappa^2 - \Delta_\Gamma)^{\alpha/2} \bar{p}(s,t) &= \bar{u}(s,t)-u_d(s,t), && \quad (s,t) \in \Gamma \times (0, T), \\ \bar{p}(s,T) &= 0, && \quad s \in \Gamma, \end{aligned} \right. \end{equation}\] with \[\begin{align} \label{zz} \tag{8} \bar{z}(s,t) = \max\left\{a(s,t),\min\left\{b(s,t),-\dfrac{1}{\mu}\bar{p}(s,t)\right\}\right\}. \end{align}\] By considering the change of variable \(\bar{q}(s,t) = \bar{p}(s,T-t)\), the fractional adjoint problem \(\eqref{eq:adjointeq}\) becomes a forward-in-time problem where the transformed adjoint state \(\bar{q}\) solves \[\begin{equation} \label{transformed_adjoint_state} \tag{9} \left\{ \begin{aligned} \partial_t \bar{q}(s,t) + (\kappa^2 - \Delta_\Gamma)^{\alpha/2} \bar{q}(s,t) &= \bar{v}(s,t)-v_d(s,t), && \quad (s,t) \in \Gamma \times (0, T), \\ \bar{q}(s,0) &= 0, && \quad s \in \Gamma, \end{aligned} \right. \end{equation}\] since \(\partial_t\bar{q}(s,t) = -\partial_t\bar{p}(s,T-t)\) and \(\bar{q}(s, 0) =\bar{p}(s,T-0)= \bar{p}(s,T)=0\). Here, \(\bar{v}(s,t) = \bar{u}(s,T-t)\) and \(v_d(s,t) = u_d(s,T-t)\).
Given \(\bar{u}\) and \(u_d\), we can time-reverse them to obtain \(\bar{v}\) and \(v_d\) and then use the same numerical scheme we use for the forward problem \(\eqref{eq:maineqoptimal}\) to solve the adjoint problem \(\eqref{transformed_adjoint_state}\). The control variable \(\bar{z}\) is then computed using \(\eqref{zz}\).
The numerical scheme for \(\eqref{eq:maineqoptimal}\) and \(\eqref{transformed_adjoint_state}\) are given by (see the Functionality page) \[\begin{align} \label{numericalscheme1} \tag{10} \bar{\boldsymbol{U}}^{k+1} = \boldsymbol{C}^{-1}\left(\sum_{k=1}^{m+1} a_k\left(\dfrac{\boldsymbol{L}\boldsymbol{C}^{-1}}{\kappa^2}-p_k\boldsymbol{I}\right)^{-1} + r\boldsymbol{I}\right) (\boldsymbol{C}\bar{\boldsymbol{U}}^k+\tau (\boldsymbol{F}^{k+1}+\bar{\boldsymbol{Z}}^{k+1})), \end{align}\] and \[\begin{align} \label{thenumericalscheme2} \tag{11} \bar{\boldsymbol{Q}}^{k+1} = \boldsymbol{C}^{-1}\left(\sum_{k=1}^{m+1} a_k\left(\dfrac{\boldsymbol{L}\boldsymbol{C}^{-1}}{\kappa^2}-p_k\boldsymbol{I}\right)^{-1} + r\boldsymbol{I}\right) (\boldsymbol{C}\bar{\boldsymbol{Q}}^k+\tau (\bar{\boldsymbol{V}}^{k+1}-\boldsymbol{V_d}^{k+1})), \end{align}\] with initial conditions \(\bar{\boldsymbol{U}}^0 = \boldsymbol{U}_0\) and \(\bar{\boldsymbol{Q}}^0 = \boldsymbol{0}\). Observe that \(\eqref{numericalscheme1}\)-\(\eqref{thenumericalscheme2}\) is a coupled problem.
For notation convenience, let
Observe that the entries of \(\boldsymbol{F}\) are given by \(\boldsymbol{F}_{ik}=(\psi^i_h, f(\cdot,t_k))_{L_2(\Gamma)}\) where \(\psi^i_h\) is the \(i\)-th basis function on a mesh with \(N_h\) nodes. By introducing a finer integration mesh with \(N_{h_{\text{ok}}}\) nodes, we can approximate \(\boldsymbol{F}_{ik}\) as \(\boldsymbol{\psi}_i^\top\boldsymbol{C}^{\text{ok}}\boldsymbol{f}_k\), where \(\boldsymbol{\psi}_i=[\psi^i_h(s_1)\dots \psi^i_h(s_{N_{h_{\text{ok}}}})]^\top\) for \(i=1\dots, N_h\), \(\boldsymbol{f}_k = [f(s_1,t_k)\dots f(s_{N_{h_{\text{ok}}}},t_k)]\top\) for \(k= 0,\dots,N\), and matrix \(\boldsymbol{C}^{\text{ok}}\) has entries \(\boldsymbol{C}^{\text{ok}}_{nm} = (\psi_{h_\text{ok}}^n,\psi_{h_\text{ok}}^m)_{L_2(\Gamma)}\) for \(n,m = 1,\dots, N_{h_\text{ok}}\). Let \(\boldsymbol{f}\) denote the matrix with columns \(\boldsymbol{f}_k\) and \(\boldsymbol{\Psi}\) denote a matrix with entries \(\boldsymbol{\Psi}_{ji}=\psi^i_h(s_j)\) for \(i=1\dots, N_h\) and \(j = 1,\dots N_{h_{\text{ok}}}\). Then \[\begin{equation} \boldsymbol{F} \approx \boldsymbol{\Psi}^\top \boldsymbol{C}^{\text{ok}} \boldsymbol{f}. \end{equation}\] Similarly, if \(\boldsymbol{U_d}\) is the matrix with entries \(\boldsymbol{U_d}^{jk} = u_d(s_j,t_k)\) for \(j = 1,\dots, N_{h_{\text{ok}}}\) and \(k = 0,\dots, N\), then the matrix \(\boldsymbol{V_d}\) can be approximated as \[\begin{equation} \boldsymbol{V_d} \approx \text{reversecolumns}(\boldsymbol{\Psi}^\top \boldsymbol{C}^{\text{ok}} \boldsymbol{U_d}). \end{equation}\]
To approximate \(\boldsymbol{F}\) and \(\boldsymbol{V_d}\), we use the values of \(f\) and \(u_d\) on the integration mesh with \(N_{h_{\text{ok}}}\) nodes. However, to approximate \(\bar{\boldsymbol{Z}}\) and \(\bar{\boldsymbol{V}}\), we do not even have access to the values of \(\bar{z}\) and \(\bar{v}\) on the coarse mesh, let alone on the fine mesh. Instead, we use their approximations on the coarse mesh, given by the matrices \(\bar{\boldsymbol{z}} = \max\{\boldsymbol{A},\min\{\boldsymbol{B},-\bar{\boldsymbol{P}}/\mu\}\}\) and \(\text{reversecolumns}(\bar{\boldsymbol{U}})\), respectively. With these at hand, we can approximate \(\bar{\boldsymbol{Z}}\) and \(\bar{\boldsymbol{V}}\) by first projecting \(\bar{\boldsymbol{z}} = \max\{\boldsymbol{A},\min\{\boldsymbol{B},-\bar{\boldsymbol{P}}/\mu\}\}\) and \(\text{reversecolumns}(\bar{\boldsymbol{U}})\) onto the fine mesh and then approximating the inner products \((\psi^i_h,\bar{z}(\cdot,t_k))_{L_2(\Gamma)}\) and \((\psi^i_h,\bar{u}(\cdot,T-t_k))_{L_2(\Gamma)}\) as \[\begin{equation} \bar{\boldsymbol{Z}} \approx \boldsymbol{\Psi}^\top \boldsymbol{C}^{\text{ok}} \boldsymbol{\Psi} \bar{\boldsymbol{z}}\quad \text{ and }\quad \bar{\boldsymbol{V}} \approx \text{reversecolumns}(\boldsymbol{\Psi}^\top \boldsymbol{C}^{\text{ok}} \boldsymbol{\Psi} \bar{\boldsymbol{U}}). \end{equation}\]
We follow the approach in Glusa and Otárola (2021) to construct an exact solution to the optimal control problem \(\eqref{eq:costfun}\)-\(\eqref{control_constraints}\). The solution to the elliptic problem \[\begin{align} \label{eq:elliptic_problem} \tag{12} (\kappa^2 - \Delta_\Gamma)^{\alpha/2} u(s)=f(s),\quad s\in\Gamma, \end{align}\] where \(u(\cdot)\) satisfies the Kirchhoff vertex conditions \(\eqref{eq:Kcond}\), is given by \[\begin{align*} u(s) = \sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}(f,e_j)_{L_2(\Gamma)}e_j(s),\quad s\in\Gamma. \end{align*}\] We choose coefficients \(\{x_j\}_{j=1}^{\infty}\) and \(\{y_j\}_{j=1}^{\infty}\) and set \[\begin{align} \label{eq:rhsfunctions} \tag{13} f(s) = \sum_{j=1}^{\infty}x_je_j(s)\quad \text{ and } \quad g(s) = \sum_{j=1}^{\infty}y_je_j(s). \end{align}\] Hence the corresponding solutions are given by \[\begin{align} \label{eq:sols} \tag{14} u(s) = \sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}x_je_j(s)\quad \text{ and } \quad v(s) = \sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}y_je_j(s). \end{align}\] We also define \(\psi(t) = \cos(t)\) and \(\phi(t) =\sin(T-t)\). Observe that \(\psi(0)=1\) and \(\phi(T)=0\). We also choose \(\mu =0.1\), \(a(s,t) = -0.5\) and \(b(s,t)=0.5\). Now we set \[\begin{align*} f(s,t) &= \psi'(t)u(s)+\psi(t)f(s) - \bar{z}(s,t) = -\sin(t)\sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}x_je_j(s)+ \cos(t)\sum_{j=1}^{\infty}x_je_j(s)- \max\{-0.5,\min\{0.5, \sin(T-t)\sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}y_je_j(s)\}\}\\ u_d(s,t) & = \psi(t)u(s)-\mu\phi'(t)v(s)+\mu\phi(t)g(s) = \cos(t)\sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}x_je_j(s)+\mu\cos(T-t)\sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}y_je_j(s)+\mu\sin(T-t)\sum_{j=1}^{\infty}y_je_j(s)\\ u_0(s) & = u(s)= \sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}x_je_j(s) \end{align*}\] The exact solution to the optimal control problem is given by \[\begin{align} \bar{u}(s,t)&=\psi(t)u(s) = \cos(t)\sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}x_je_j(s)\\ \bar{p}(s,t)&=-\mu\phi(t)v(s)= - \mu\sin(T-t)\sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}y_je_j(s)\\ \bar{z}(s,t)&=\max\{a(s,t),\min\{b(s,t), \phi(t)v(s)\}\}= \max\{-0.5,\min\{0.5, \sin(T-t)\sum_{j=1}^{\infty}\lambda_j^{-\alpha/2}y_je_j(s)\}\} \end{align}\]
Observe that the choice of \(f\) decouples the systems \(\eqref{numericalscheme1}\)-\(\eqref{thenumericalscheme2}\).
At \(t=0\), we get \[\begin{align*} \bar{u}(s,0) = \psi(0)u(s) = u(s) = u_0(s), \end{align*}\] while for \((s,t) \in \Gamma \times (0, T)\), we obtain \[\begin{align*} \partial_t \bar{u}(s,t) +(\kappa^2 - \Delta_\Gamma)^{\alpha/2} \bar{u}(s,t) = \psi'(t)u(s) + \psi(t) (\kappa^2 - \Delta_\Gamma)^{\alpha/2} u(s) = \psi'(t)u(s) + \psi(t)f(s) = f(s,t) + \bar{z}(s,t). \end{align*}\] Similarly, at \(t=T\), we get \[\begin{align*} \bar{p}(s,T) = -\mu\phi(T)v(s) = 0, \end{align*}\] while for \((s,t) \in \Gamma \times (0, T)\), we obtain \[\begin{align*} -\partial_t \bar{p}(s,t) +(\kappa^2 - \Delta_\Gamma)^{\alpha/2} \bar{p}(s,t) = \mu\phi'(t)v(s) -\mu \phi(t) (\kappa^2 - \Delta_\Gamma)^{\alpha/2} v(s) = \mu\phi'(t)v(s) -\mu \phi(t)g(s) \end{align*}\] and \[\begin{align*} \bar{u}(s,t) - u_d(s,t) = \psi(t)u(s) - \psi(t)u(s) + \mu\phi'(t)v(s) - \mu\phi(t)g(s) = \mu\phi'(t)v(s) -\mu \phi(t)g(s). \end{align*}\] This shows that \(\bar{u}\) and \(\bar{p}\) satisfy \(\eqref{eq:maineqoptimal}\) and \(\eqref{eq:adjointeq}\).
In this section, we implement the numerical scheme for the optimal control problem \(\eqref{eq:costfun}\)-\(\eqref{control_constraints}\) on the tadpole graph. We start by defining the parameters and constructing the exact solution. We them use \(\eqref{numericalscheme1}\)-\(\eqref{thenumericalscheme2}\) to solve the forward and adjoint problems. Finally, we compute the control variable \(\bar{z}\) using \(\eqref{zz}\).
# Parameters
T_final <- 2
time_step <- 0.02
h <- 0.02
kappa <- 15
alpha <- 1.8
m = 4
beta <- alpha/2
mu <- 0.1
a <- - 0.5*100
b <- 0.5*100
N_finite = 4 # choose even
adjusted_N_finite <- N_finite + N_finite/2 + 1
# Coefficients for f and g
coeff_elliptic_g <- rep(0, adjusted_N_finite)
coeff_elliptic_g[7] <- 100*100
coeff_elliptic_f <- rep(0, adjusted_N_finite)
coeff_elliptic_f[7] <- -10*100
# time and spatial mesh
time_seq <- seq(0, T_final, length.out = ((T_final - 0) / time_step + 1))
graph <- gets.graph.tadpole(h = h)
graph$compute_fem()
G <- graph$mesh$G
C <- graph$mesh$C
L <- kappa^2*C + G
I <- Matrix::Diagonal(nrow(C))
# Compute the eigenvalues and eigenfunctions
eigen_params <- gets.eigen.params(N_finite = N_finite,
kappa = kappa,
alpha = alpha,
graph = graph)
EIGENVAL_MINUS_ALPHA <- eigen_params$EIGENVAL_MINUS_ALPHA
EIGENFUN <- eigen_params$EIGENFUN
# Construct the right hand side functions f and g for the elliptic problem
elliptic_f <- as.vector(EIGENFUN %*% coeff_elliptic_f)
elliptic_g <- as.vector(EIGENFUN %*% coeff_elliptic_g)
# Construct the corresponding elliptic solution u and v
elliptic_u <- as.vector(EIGENFUN %*% (coeff_elliptic_f * EIGENVAL_MINUS_ALPHA))
elliptic_v <- as.vector(EIGENFUN %*% (coeff_elliptic_g * EIGENVAL_MINUS_ALPHA))
A <- matrix(a, nrow = nrow(C), ncol = length(time_seq))
B <- matrix(b, nrow = nrow(C), ncol = length(time_seq))
# Construct the optimal variables
psi <- cos(time_seq)
psi_prime <- - sin(time_seq)
phi <- sin(T_final - time_seq)
phi_prime <- - cos(T_final - time_seq)
u_bar <- elliptic_u %*% t(psi)
p_bar <- - mu * (elliptic_v %*% t(phi))
z_bar <- pmax(A, pmin(B, - p_bar / mu))
# Construct the data for the problem
f <- elliptic_u %*% t(psi_prime) +
elliptic_f %*% t(psi) -
z_bar
u_d <- elliptic_u %*% t(psi) -
mu * (elliptic_v %*% t(phi_prime)) +
mu * (elliptic_g %*% t(phi))
# Construct the fractional operator, which is shared for the forward and adjoint problems
my_op_frac <- my.fractional.operators.frac(L, beta, C, scale.factor = kappa^2, m = m, time_step)
# Construct the integration mesh
overkill_graph <- gets.graph.tadpole(h = 0.0001)
overkill_graph$compute_fem()
overkill_weights <- overkill_graph$mesh$weights
overkill_EIGENFUN <- gets.eigen.params(N_finite = N_finite,
kappa = kappa,
alpha = alpha,
graph = overkill_graph)$EIGENFUN
# Construct the right hand side functions f and g for the elliptic problem on the integration mesh
overkill_elliptic_f <- as.vector(overkill_EIGENFUN %*% coeff_elliptic_f)
overkill_elliptic_g <- as.vector(overkill_EIGENFUN %*% coeff_elliptic_g)
# Construct the corresponding elliptic solution u and v on the integration mesh
overkill_elliptic_u <- as.vector(overkill_EIGENFUN %*% (coeff_elliptic_f * EIGENVAL_MINUS_ALPHA))
overkill_elliptic_v <- as.vector(overkill_EIGENFUN %*% (coeff_elliptic_g * EIGENVAL_MINUS_ALPHA))
# Construct the projection matrix
Psi <- graph$fem_basis(overkill_graph$get_mesh_locations())
R <- t(Psi) %*% overkill_graph$mesh$C
# overkill_A <- matrix(a, nrow = length(overkill_weights), ncol = length(time_seq))
# overkill_B <- matrix(b, nrow = length(overkill_weights), ncol = length(time_seq))
#
#
# z_bar_overkill <- pmax(overkill_A, pmin(overkill_B, overkill_elliptic_v %*% t(phi)))
# f <- overkill_elliptic_u %*% t(psi_prime) +
# overkill_elliptic_f %*% t(psi) -
# z_bar_overkill
# U_d <- overkill_elliptic_u %*% t(psi) -
# mu * (overkill_elliptic_v %*% t(phi_prime)) +
# mu * (overkill_elliptic_g %*% t(phi))
# u_0 <- elliptic_u
f <- (overkill_EIGENFUN[, 3] + overkill_EIGENFUN[, 5]) %*% t(phi)
U_d <- (overkill_EIGENFUN[, 2] + overkill_EIGENFUN[, 4] + + overkill_EIGENFUN[, 6]) %*% t(psi)
u_0 <- EIGENFUN[, 7]
F_proj <- R %*% f
V_d <- reversecolumns(R %*% U_d)tol <- 1e-12
maxit <- 25
verbose <- TRUE
weights <- graph$mesh$weights
res <- solve_coupled_system_multi_tol(
my_op_frac = my_op_frac,
time_step = time_step,
time_seq = time_seq,
u_0 = u_0,
F_proj = F_proj,
V_d = V_d,
Psi = Psi,
R = R,
A = A,
B = B,
mu = mu,
weights = weights,
tol = tol,
maxit = maxit,
verbose = verbose,
true_sol = list(
u_bar = u_bar,
p_bar = p_bar,
z_bar = z_bar
))## iter 1: rel(Z) = 2.720e+02, rel(U) = 1.000e+00, rel(P) = 1.000e+00## iter 2: rel(Z) = 1.405e-01, rel(U) = 1.501e+00, rel(P) = 1.405e-01## iter 3: rel(Z) = 7.276e-05, rel(U) = 7.752e-04, rel(P) = 7.276e-05## iter 4: rel(Z) = 3.784e-08, rel(U) = 4.023e-07, rel(P) = 3.784e-08## iter 5: rel(Z) = 1.977e-11, rel(U) = 2.097e-10, rel(P) = 1.977e-11## iter 6: rel(Z) = 1.037e-14, rel(U) = 1.097e-13, rel(P) = 1.037e-14## [1] TRUE## [1] 6